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1  Overview 

This  report  summarizes  results  obtained  within  the  last  Funding  Period  of  the  Project  “Small- 
Molecule  Reactions  Relevant  to  the  Hypersonic  Flight  Regime”.  Part  of  the  work  has  been 
published[l,  2]  and  other  work  is  currently  prepared  for  publication. [3] 

2  Equilibrium  Rate  Coefficients  from  Atomistic  Simulations: 
The  0(3P)  +  NO(2n)  ^  02{X^^g)  +  N{^S)  Reaction  at 
Temperatures  Relevant  to  the  Hypersonic  Flight  Regime 

Reactions  involving  nitrogen-  and  oxygen-containing  small  molecules  occur  in  a  wide  range 
of  processes.  As  an  example,  NO2  -  which  can  decay  into  NO  -|-  O  or  O2  -|-  N  -  plays  a  major 
role  in  atmospheric  chemistry,  as  a  smog  constituent  and  in  combustion  processes.  [4,  5,  6]  The 
reaction  characteristics  changes  from  formation  of  NO2  to  exchange  and  dissociation  reactions 
at  higher  temperatures  (1000-20000  K).  Above  20000  K,  it  is  expected  that  the  reactions  will 
be  dominated  by  the  complete  dissociation  to  the  atomic  species.  In  the  hypersonic  flight 
regime  -  which  is  the  physical  situation  of  interest  in  the  present  work  -  the  chemistry  will  be 
in  the  intermediate  case.  Near  the  surface  of  such  vehicles  highly  non-equilibrium  conditions 
will  be  present  with  vibrational  and  rotational  temperatures  independently  reaching  several 
thousand  Kelvin.  [7]  The  gas-phase,  surface  reactions  and  energy  transfer  at  these  tempera¬ 
tures  are  essentially  uncharacterized  and  the  experimental  methodologies  capable  of  probing 
them  are  not  well  established.  In  this  respect,  theoretical  and  validated  computational  ap¬ 
proaches  become  a  valuable  complement  to  experiments. 

In  order  to  model  chemical  kinetics  under  such  non-equilibrium  conditions,  state-  and  temperature- 
dependent  rate  coefficients  for  individual  reactions  are  required.  However,  direct  experimen¬ 
tation  under  such  conditions  is  usually  not  possible.  In  the  present  work,  both  forward  and 
reverse  reaction  rate  coefficients  for  0(^P)  -|-  NO(^n)  o  02{X^T,~)  +  N(^S)  are  calculated 
from  molecular  dynamics  (MD)  simulations  for  moderate  to  high  temperatures  (1000  to  20000 
K)  and  compared  with  results  from  experimentally  derived  thermodynamics  quantities  from 
the  NASA  CEA  (NASA  Chemical  Equilibrium  with  Applications)  database.  [8,  9]  The  value 
of  such  a  validated  computational  approach  is  that  it  can  be  applied  generically  to  a  wide 
range  of  bimolecular  reactions  whenever  direct  experiments  are  not  available,  with  consis¬ 
tency  of  the  computations  depending  on  a  comparison  of  independent  dynamic  and  energetic 
approaches. 

A  schematic  representation  of  the  relevant  species  for  the  0(^P)  -|-  NO(^n)  <-)•  02(^E“)  -|- 
N(^P)  reaction  is  given  in  Figure  I.  NO2  formation  has  a  barrier  of  1.3  eV  from  the  N-I-OI02 
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side,  whereas  from  the  NOl+02  side  it  is  a  barrierless  process.  In  order  to  reach  the  N+0102 
side  from  the  NOl+02  side,  2.6  eV  are  required  including  the  barrier.  An  additional  5.6  eV 
are  required  for  dissociation  of  O2  and  6.9  eV  for  the  dissociation  of  NO  to  the  atomization 
state  of  N+01+02. 


O  +  O  +  N 


r23=  4.12  a„ 

01  N 


Figure  1:  Relevant  stages  in  the  NOl+02  reaction.  The  energies  were  derived  from  the  ab 
initio  calculations  in  Ref.  [10]. 

Even  when  the  excited  states  can  approach  the  ground  state  PES  for  certain  configurations 
(Eigure  2),  their  influence  can  be  considered  to  be  minor.  The  first  excited  state  is  dissocia¬ 
tive  and  about  6  eV  are  required  to  access  the  second  and  third  excited  states  (see  Figure 
2).  For  the  Maxwell-Boltzmann  distribution  at  20000  K,  translational  energies  beyond  6  eV 
are  accessed  with  little  probability  (0.06).  This  provides  an  upper  bound  since  a  thermal  dis¬ 
tribution  requires  additional  energy  to  be  distributed  into  rotational  and  vibrational  motion. 
Furthermore,  lifetimes  of  the  system  at  the  second  and  third  excited  states  are  presumably 
small  since  the  wells  are  shallow.  This  causes  the  system  to  rapidly  return  to  the  ground 
state.  [11]  The  expected  effect  of  including  electronically  excited  states  is  to  somewhat  lower 
the  rate  coefficients  due  to  the  topology  of  the  excited  PES  (see  Figure  2  which  exhibits 
shallow  minima  that  potentially  stabilize  intermediates  for  short  times).  Therefore,  in  this 
work  we  first  consider  only  the  ground  state  PES  since  it  dominates  the  reaction  dynamics 
under  the  relevant  physicochemical  conditions. 

The  equilibrium  constant  as  a  function  of  temperature  for  a  chemical  reaction  at  equilibrium 
is  defined  as,  [12] 

K(T)  =  (1) 

where  kj^{T)  and  k-{T)  are  the  rate  coefficients  for  the  forward  and  reverse  reactions,  respec- 
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N-0  distance  (a  J 


Figure  2:  MRCI+Q/cc-pVQZ  energies  (symbols)  for  the  O-NO  coordinate  with  the  other 
N-0  distance  fixed  at  r^o  =  2.25  oq  for  linear  NO2  arrangement.  The  solid  lines  are  RKHS 
interpolants  for  each  electronic  state.  The  energy  gap  between  the  ground  and  third  electronic 
state  has  been  indicated  by  a  vertical  arrow.  The  panel  on  the  right  hand  side  reports  the 
probability  distribution  function  P{E)  for  the  total  energy  E. 

tively.  The  rate  coefficient  for  the  forward  reaction  (A:+(T))  has  been  recently  determined 
from  atomistic  simulations.  [10]  The  same  methodology  is  used  here  for  the  calculation  of  the 
rate  coefficient  for  the  reverse  reaction  k-{T).  Following  the  previous  work,  three  steps  are 
involved:  (1)  the  construction  of  the  PES  for  the  ground  state  of  the  NO2  molecule  based 
on  high-level  ab  initio  calculations  and  its  representation  with  a  reproducing  kernel  Hilbert 
space  (RKHS)  method[13]  combined  with  Legendre  polynomials;  (2)  quasi  classical  trajectory 
(QCT)  calculations  to  study  the  adiabatic  reaction  dynamics,  and  (3)  calculation  of  the  rate 
coefficients  for  the  different  exit  channels  using  an  Importance  Sampling  Monte  Carlo  method. 

The  thermal  rate  coefficient  from  trajectory  calculations  is  determined  from[14] 

I -  00 

k{Tt,T„,  Te)  =  J ^  [ a{E,- r„)Ece-^‘®MEc,  (2) 

g[%)  V  TTg  J 

*  0 

where  Pt  =  and  ks  is  the  Boltzmann  constant.  Ft,  Trv,  Te  are  the  translational  tem¬ 
perature  of  NO  and  O,  the  rovibrational  temperature  of  NO  and  the  electronic  temperature 
of  NO,  respectively,  g{Te)  is  the  electronic  degeneracy  factor, [11,  15]  g  is  the  reduced  mass 
of  NOl  and  O2,  respectively,  and  a{Ec]T„)  is  the  integral  cross  section  as  a  function  of  the 
collision  energy,  E^  and  In  our  studies  it  was  assumed  that  Tt  =  Trv  =  Te.  In  other  words, 
the  various  degrees  of  freedom  in  the  reactants  are  in  thermal  equilibrium.  If  no  subscript 
is  shown  for  T  then  it  corresponds  to  the  common  temperature.  In  a  similar  manner,  the 
reaction  cross  section  a{Ec‘,T„)  for  a  given  collision  energy  E^  is 
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(t{Ec,T„) 


■i^max  3max{'^) 

t?=0  j=0 

fmax  ima3:('^^) 

yy  (2j  +  i)e“/^rvE„j 

■!;=0  j=0 


(3) 


where  cj„j(£'c)  is  the  (v,j)— state  dependent  cross  section  at  collision  energy  E^.  The  energy 
E^j  of  a  rovibrational  state  {v,j)  is  calculated  according  to  a  Morse  oscillator  model  for  the 
NO  molecule. [16]  The  cross  section  as  an  integral  of  the  opacity  function  Pvj{b',  E^)  for  given 
Ec  and  rovibrational  state  (u,j)  is 


ayj{Ec)  =  J  Pvj{b;Ec)27rbdb  (4) 

0 

The  integral  in  Eq.  2  can  be  calculated  by  using  an  Importance  Sampling  Monte  Carlo 
scheme. [17]  For  this,  the  vibrational  and  rotational  quantum  numbers  v  and  j,  and  the  colli¬ 
sion  energy  (Ec)  are  sampled  from  the  Boltzmann- weighted  probability  distribution  function 


and 


respectively. 


Pvj  (Trv) 


(2j  -I-  l)e 

Vmax  jrnaxiy') 

yy  yy  (2/ 

v'=0  j'=0 


piE,)dE,  =  (3^E,e-^^^^dE„ 


(5) 


(6) 


For  the  determination  of  the  upper  limit  (6max)  of  the  relevant  sampling  interval  for  the  im¬ 
pact  parameter,  b,  preliminary  tests  were  carried  out  at  low  temperatures  (100  and  200  K). 
As  NO2  was  not  formed  for  6  >  26  ao,  the  impact  parameter  was  therefore  uniformly  sampled 
from  0  up  to  6max  =  26  ao.  As  the  opacity  curves  in  Figure  3  show,  this  upper  limit  is  very 
conservative  and  ensures  convergence  of  the  results  with  respect  to  b.  However,  large  b  is 
required  for  low-energy  collisions.  [10] 


After  evaluating  the  integral  in  Eq.  2  over  v,  j  and  E^  the  resulting  expression  for  k(T)  is 


k{T) 


8  2Trb 


-^reac 

'max  1 

>  bi, 


TT/U/?  g(T)Ntot  ^ 


(7) 


where  Aj-eac  and  Atot  are  the  number  of  reactive  and  the  total  number  of  trajectories,  respec¬ 
tively,  and  bi  is  the  impact  parameter  of  reactive  trajectory  i.  The  convergence  of  the  integral 
in  Eq.  2  was  monitored  by  the  decrease  of  the  Monte  Carlo  error.  [10]  Choosing  eigenenergies 
from  a  Morse-oscillator  for  the  NO  molecule  was  done  for  convenience.  To  quantify  the  ex¬ 
pected  difference  in  using  a  kernel-interpolated  1-dimensional  potential  for  the  NO  oscillator, 
the  bound  states  from  the  Morse-fit  and  the  RKHS  were  determined  up  to  u  =  6  (correspond¬ 
ing  to  ~  12000  cm“^).  The  correlation  between  the  two  sets  of  eigenenergies  is  better  than 


DISTRIBUTION  A.  Approved  for  public  release:  distribution  unlimited. 


Figure  3:  Thermal  opacity  functions  for  1000  and  10000  K  on  the  upper  and  lower  panels,  re¬ 
spectively.  The  calculated  points  for  N02  formation  in  the  high  pressure  limit  (open  squares), 
oxygen  exchange  channel  (triangles)  and  the  ”0102-|-N”  channel  (open  diamonds)  have  been 
represented.  Error  bars  are  indicated  in  black. 


0.999  with  a  slope  of  1.0003.  Hence,  over  the  relevant  energy  range  the  initial  conditions  from 
a  Morse-oscillator  model  are  expected  to  closely  reflect  those  from  using  the  more  accurate 
kernel-interpolated  potential. 

The  rate  coefficients  for  the  reverse  {k-{T))  and  forward  reaction  (A:+(T))  derived  directly 
from  trajectory  calculations [10]  in  the  present  work  are  reported  in  Figure  4  for  temperatures 
between  300  and  20000  K.  An  inset  with  results  for  the  range  from  1000  to  5000  K  has  also 
been  included.  The  curve  from  this  work  represents  a  fit  to  a  modified  Arrhenius  equation 
(see  Eq.  8). [18]  A  total  of  10000  trajectories  was  run  for  the  NO-l-0  reaction  whereas  only 
5000  were  necessary  for  the  reverse  reaction  to  obtain  a  rate  coefficient  with  a  relative  error 
of  better  than  10  %.  It  is  found  that  a  directly  measured  data  point  at  3000  K  (green  star) 
agrees  very  favorably  with  the  atomistic  simulations  (magenta)  for  /c_|_(r).[10] 

k{T)  =  (8) 

where  A  is  the  pre-exponential  factor  and  Ea  is  the  activation  energy.  The  parameters  result¬ 
ing  from  the  fit  are  A  =  9.35  x  10®  cm®s“^mol“^K“”,  Ea  =  1.88  x  10^  kcal/mol  and  n  =  0.93. 

Forward  rate  coefficients  {k+{T))  have  been  reported  in  previous  works  most  of  which  were 
based  on  computation.  Valli  et  al  reported  A:+(T)  from  quasi  classical  trajectory  calculations 
using  the  ground  state  PES  of  the  NO2  molecule  for  temperatures  ranging  from  300  to  500 
K.[19]  This  data  has  been  extended  to  higher  temperatures  by  means  of  a  fit  to  Eq.  8.  The 
resulting  parameters  are  A  =  1.65  x  10^®  cm®s“^mol“^K“”,  Ea  =  2.13  x  10^  kcal/mol  and 
n  =  0.  In  the  review  of  reactions  relevant  to  Combustion  Chemistry  by  Baulch  et  al.  A:+(T) 
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for  temperatures  ranging  from  1000  to  5000  K  are  reported  and  fit  to  Eq.  8  with  parameters 
A  =  6.87  X  10®  cm®s“^mol“^K“"',  Ea  =  1.9  x  10^  kcal/mol  and  n  =  1.13. [20] 

Bose  et  al  reported  rate  coefficients  for  the  reverse  reaction  (A:_(T))  based  on  QCT  on  the  “^A' 
and  '^A!  PESs  for  temperatures  ranging  from  1000-14000  K.[21]  In  order  to  obtain  the  forward 
reaction  rate  coefficients  using  Eq.l,  these  data  were  combined  with  CEA  data[8,  9]  for  the 
equilibrium  rates  {K(T))  (available  from  300-20000  K).  These  results  were  extended  over  the 
whole  range  of  temperatures  from  300-20000  K  by  fitting  Eq.8  which  yielded  the  following 
parameters  A  =  2.49  x  10®  cm®s“^mol“^K“”',  Ea  =  0.4  x  10^  kcal/mol  and  n  =  1.18. 

At  lower  temperatures  (below  1000  K)  experimental  data  is  available  for  which  compares 
very  favourably  with  the  computed  values  as  well.  At  sufficiently  low  temperatures  quan¬ 
tum  effects  (effect  of  zero  point  energy  and  tunneling)  will  become  important.  A  recent 
study  of  the  proton  transfer  reaction  in  malonaldehyde  has  established  that  zero  point  effects 
can  play  an  important  role  whereas  tunneling  probably  only  becomes  important  at  the  low¬ 
est  temperatures.  [22]  However,  this  temperature  range  is  not  of  primary  importance  for  the 
present  work. 


0  5000  10000  15000  20000 

T,K 


Figure  4:  Rate  coefficients  k-{T)  (blue  open  triangles)  and  k^{T)  (red  open  triangles)  for 
temperatures  between  300  and  20000  K.  An  inset  with  results  for  the  range  from  1000  to  5000 
K  has  been  included.  In  the  legend  [a]=Ref.[19],  [b]=Ref.[20],  [c]=Ref.[21],  [d]=Ref.[10]  and 
[e]=Ref.[23]. 

The  equilibrium  rates  from  the  present  work  obtained  directly  from  the  trajectory  calcula¬ 
tions  are  compared  to  the  equilibrium  constants  calculated  from  the  CEA  database  [8,  9,  24] 
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in  Figure  5  for  temperatures  between  5000  and  20000  K.  For  temperatures  below  6000  K,  the 
CEA  functions  are  fit  to  available  experimental  data  whereas  above  6000  K,  the  CEA  results 
have  been  obtained  from  usual  statistical  thermodynamic  expressions  for  equilibrium  rates, 
based  on  spectroscopic  constants  of  the  atoms  and  molecules  involved.  The  results  in  Eigure  5 
show  that  CEA  equilibrium  constants  are  somewhat  underestimated  for  temperatures  below 
10000  K  and  agree  almost  quantitatively  for  higher  temperatures. 


Eigure  5:  Equilibrium  rate  coefficients  K{T)  for  temperatures  between  5000  and  20000  K. 
Results  from  CEA  are  also  included  for  comparison.  [8,  9] 

In  conclusion,  the  equilibrium  rate  coefficient  for  the  0(^P)  +  NO(^n)  -f-)-  02{X^T,~)  +  N(^S) 
reaction  was  calculated  directly  from  molecular  dynamics  simulations  on  accurate  PESs  fitted 
to  high-level  electronic  structure  calculations.  The  K{T)  thus  obtained  compare  favorably 
with  experimental  measurements  from  CEA  over  a  wide  temperature  range  extending  up  to 
20000  K.[8,  9,  24]  The  influence  of  higher  electronic  states  for  higher  temperature  is  minor. 
The  forward  rate  coefficients  from  this  work  have  been  validated  against  theoretical  results 
from  previous  works.  [19,  20,  21]  and  the  present  ansatz  provides  a  generic  computational 
framework  for  obtaining  such  data  from  computation. 

3  Reaction  Dynamics  for  the  N(^S)+NO(^n)  System  at  Tem¬ 
peratures  Relevant  to  the  Hypersonic  Flight  Regime 

Nictric  oxide  (NO)  plays  a  major  role  in  the  chemistry  near  the  surface  of  vehicles  for  at¬ 
mospheric  re-entry  in  the  hypersonic  flight  regime  [25].  Therefore,  reactions  of  nitric  oxide 
with  the  most  common  atmospheric  gases,  e.g.  N  and  O,  must  be  considered  in  the  design 
of  spacecraft.  Modelling  the  thermochemical  phenomenon  in  the  re-entry  flows  requires  the 
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knowledge  of  rate  coefficients  of  such  reactions  at  very  high  temperatures.  At  such  extreme 
conditions  quantitative  experiments  are  difficult  and  expensive.  Consequently,  computational 
studies  are  vital  for  characterizing  the  reactivity  or  energy  transfer  in  reactive  or  nonreactive 
encounters  of  such  systems. 

Recently,  a  study  of  NO  formation  in  hypersonic  boundary  layers  [25]  around  space  vehicles 
during  re-entry  showed  that  the  inclusion  of  the  NO  data  in  the  models  modified  the  mass 
fraction  and  the  heat  flux  substantially  in  the  boundary  layers.  The  ground  state  of  N2O 
has  been  studied  intensely  both,  theoretically  and  experimentally.  [26,  27,  28,  29,  30,  31,  32, 
33,  34]  However,  the  potential  energy  surfaces  (PESs)  that  connect  NO(X^n)-|-N('^5)  with 
N2(A^S)-|-0(^P)  are  in  the  triplet  manifold,  see  Figure  6.  Several  experimental  studies  of 
the  NO(X^n)-|-N(‘^5)— >N2(A^S)-|-0(^P)  reaction  have  been  performed  over  the  past  decades 
[35,  36,  34,  37,  38,  39]  Nevertheless,  the  agreement  between  the  experimental  results  is  poor, 
e.g.  at  100  K  Fox[40]  derived  a  rate  of  3.60x10“^^  cm^molecule“^s“^  while  Bergeat  et  al  [41] 
report  a  value  of  4.11x10“^^  cm^molecule“^s“^. 


R  (A)  R  (A) 


Figure  6:  Schematic  representation  of  the  Coou  and  Cg  {NON— angle  near  130°)  partial  adi¬ 
abatic  correlation  diagrams  for  N2O.  Data  taken  from  Ref.  [42].  The  Boltzmann  distribution 
at  1000  K  (red  line)  and  10000  K  (blue  line)  for  the  initial  NO(X^n)  are  represented  together 
with  the  lowest  vibrational  states. 


Theoretical  studies  of  the  scattering  on  the  adiabatic  electronic  states  have  been  done  on  the 
and  ^A"  surfaces.  [43,  44,  45,  46]  The  PESs  of  Garmallo  et  al  [45]  have  been  considered 
up  recently  as  the  best  reported  in  the  literature  [25].  The  analytical  representation  of  these 
surfaces  is  based  on  a  many-body  expansion  derived  from  a  grid  of  ab-initio  data  computed 
with  the  CASPT2  method.  Very  recently,  new  PESs  for  the  two  lowest  triplet  states  were 
published.  [47]  These  surfaces  were  calculated  at  the  multi  reference  configuration  interaction 
(MRCI)  level  and  showed  remarkable  differences  compared  to  the  PESs  of  Garmallo  et  al  [45]. 
The  global  minimum  of  the  ^A"  was  found  to  be  0.23  eV  deeper  than  in  the  previous  work. 
For  the  ^A'  state  several  new  minima  and  transition  states  were  found.  Therefore,  MRGI 
calculations  are  required  for  computing  meaningful  rates. 

The  new  available  PESs[47]  are  limited  to  study  the  N2-I-O — )•  NO-I-N  forward  reaction.  This 
is  due  to  the  R~^  dependence  in  the  long  range  (dispersion  and  induction  interaction)  is  not 
included  explicitly  in  the  fit.  The  N2-I-O  reaction  is  highly  endothermic,  and  the  dynamics  in 
this  direction  is  very  expensive  computationally.  Very  recently,  the  same  group  presented  a 

DISTRIBUTION  A.  Approved  for  public  release:  distribution  unlimited. 

8 


quasiclassical  trajectory  study  of  the  N2+O — t-NO+N  using  their  surfaces  [48].  In  that  work, 
they  concluded  that  the  reaction  occurs  predominantly  in  the  surface.  However,  no  rate 
coefficients  necessary  in  the  models  of  the  atmospheric  re-entry  were  determined. 

In  the  present  work  new  PESs  for  the  and  ^A''  states  of  N2O  are  presented.  These  PESs 
are  used  to  study  the  NO(A^n)-|-N(^S')  collision  and  computed  the  thermal  rate  coefficients 
at  high  temperatures.  In  the  next  section,  the  methods  used  in  the  development  of  the  PESs 
are  described,  and  also  a  brief  description  scattering  calculations  is  included.  The  results  are 
discussed  in  Section  III,  while  the  main  conclusions  are  in  Section  IV. 

3.1  Methods 
Ab-initio  calculations 

The  ab-initio  potential  energies  of  the  triplet  states  of  N2O  were  computed  at  the  multirefer¬ 
ence  configuration  interaction  (MRCI)  level  including  the  Davidson  correction  (MRCI-|-Q),[49] 
preceded  by  complete  active  space  self-consistent  field  (CASSCE)  calculations.  The  standard 
Dunning’s  correlation  consistent  polarized  triplet-zeta  (cc-pVTZ)  basis  set  [50]  were  used. 
Also,  a  set  of  points  was  computed  using  the  coupled  cluster  method  with  single  and  double 
electronic  excitations  and  perturbative  treatment  of  triple  excitation  (CCSD(T))  for  ascertain¬ 
ing  the  quality  of  the  MRCI-I-Q  calculations  in  the  long  range  part  {R  — )■  00).  All  electronic 
structure  calculations  were  performed  with  the  Molpro  package  [51]. 

In  both  electronic  states,  ^A'  and  ^A",  the  asymptotic  N2-I-O  and  NO-I-N  channels  were 
treated  separately.  In  each  case,  energies  on  a  three-dimensional  grid  for  the  respective  Jacobi 
coordinates  {R,  r  and  6)  were  computed.  Eor  N2-I-O,  the  angles  vary  between  0°  and  90°.  The 
angular  interval  was  extended  up  to  180°,  taking  into  account  the  symmetry  of  the  system. 
In  this  way,  the  number  of  angular  points  is  the  same  that  for  NO-I-N,  which  is  necessary  for 
the  fitting.  In  total,  11  angular  points  were  taken  into  account  in  each  channel.  MRCI-I-Q 
calculations  for  the  ^A'  state  of  the  N2-I-O  asymptote  were  carried  out  for  19  distances  along 
R  {R  ^  [1.1,10.0]  A),  and  13  points  along  r  (r  ^  [0.90,1.47]  A)  whereas  for  the  NO-I-N 
asymptote  they  included  R  £  [1.0, 10.0]  A  and  r  £  [0.90, 1.55]  A.  Eor  the  ^A"  state  the  same 
R  intervals  were  used  whereas  for  the  N2-I-O  channel  r  G  [0.90, 1.55]  A  and  r  £  [0.90, 1.58]  A 
for  the  NO-I-N  channel. 


Representation  of  the  PESs 

All  PESs  were  represented  by  using  a  reproducing  kernel  Hilbert  space  approach  (RKHS).[52] 
The  potential  energy  for  each  channel  is 

M 

l/(x)  =  afc<5(xfc,x)  (9) 

k=l 

where  M  is  the  number  of  ab-initio  points,  x  is  the  vector  of  the  internal  coordinates 
X  =  (R,r,z),  with  2;  =  corresponds  to  a  grid  point  Xk  =  {Rk,rk,  Zk)- 

Q(xk,x)  is  the  reproducing  kernel  for  a  multidimensional  RES,  which  corresponds  to  the 
multiplication  of  three  1-dimensional  reproducing  kernels 


(5(xk,x)  =  q‘^'^{Rk,  R)q^'^{rk,r)q'^{zk,  z) 


(10) 
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The  coefficients  ak  are  determined  from  solving  the  linear  system 


/  (5(xi,xi)  Q(xi,X2)  •••  (5(xi,xm)  \ 

/  ai  \ 

(5(x2,Xi)  Q(x2,X2)  •••  (5(x2,Xm) 

02 

\  (5(xm,xi)  Q(xm,X2)  •••  Q(xm,xm)  / 

\  O'  M  / 

(  ^(xi)  \ 

^(X2) 


V  ^(xm)  / 


(11) 


By  construction,  reproducing  kernels  decay  to  zero  for  i?  — )■  oo.  Hence,  for  having  a  correct 
representation  of  the  asymptotic  behaviour  of  the  kernel  interpolation,  the  energy  at  large 
atom-diatom  distance  — )•  oo,r,6)  was  subtracted  from  the  values  of  the  grid,  such  that 

each  cut  r,  6)  dissociates  to  zero.  The  potential  r,  6)  is  then  represented  as  an 

RKHS.  The  asymptotic  energies  — )•  oo,r,6),  considered  isotropic,  are  also  fitted  using 

a  1-dimensional  kernel.  The  total  energy  is  the  sum  of  the  energies  of  both  interpolations 
V {R,  r,  6)  =  r,9)  +  V{R  ^  oo,  r,  9).  The  procedure  is  similar  to  that  used  for  Ar-|-N^ 

[2]. 


The  reactive  PES  was  constructed  by  smoothly  connecting  the  three  PESs  (AB-I-C,  AC-I-B, 
and  BC-I-A),  each  represented  as  a  RKHS  (Eq.  9)  with  a  distance-dependent  switching 
function,  as  was  previously  done  for  N02.[53,  1]  The  global  PES  is  therefore 

3 

^(y)  =  (12) 

i=l 

where  the  y  is  the  vector  of  the  three  inter-atomic  distances  of  the  N2O  molecule,  and  Vi  are 
the  potential  energies  of  the  three  asymptotic  PESs.  The  switching  function  is 

p-iVi/a)* 

WiiVi)  =  -  (13) 

i=i 

where  y*  are  either  of  the  inter-atomic  distance  and  the  parameter  a  was  set  to  1.0  for  the  ^A' 
PES  and  1.2  in  the  ^A''  surface.  This  parameter  determines  the  range  over  which  the  PESs 
are  mixed. 


Furthermore,  an  independent  grid  of  points  was  computed  for  describing  the  diatomic  molecules 
N2  and  NO,  using  the  same  method  that  in  the  calculation  of  the  global  PES.  These  points 
were  fitted  to  the  Morse  potential  function  using  a  least  squares  procedure.  Table  1  summa¬ 
rizes  the  properties  of  both  molecules.  The  dissociation  energy  is  lower  than  the  experimental 
data,  however,  this  should  not  affect  the  calculations  as  the  full  dissociation  of  N2O  into 
N-|-N-|-0  is  not  considered  here.  These  Morse  parameters  are  then  used  for  generating  initial 
conditions  for  the  quasiclassical  trajectories  (QCT)  studies. 
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Table  1;  Morse  parameters  and  diatomic  data  of  N2  and  NO. 


Te 

(A) 

De 

(eV) 

/3 

(A-i) 

We 

(cm“^) 

UJ^Xq 

(cm“^) 

N2 

calc. 

1.105 

9.385 

2.799 

2389.579 

18.858 

expt. 

1.098 

9.905 

2358.570 

14.324 

NO 

calc. 

1.158 

6.234 

2.942 

1982.032 

19.532 

expt. 

1.151 

6.614 

1904.200 

14.075 

Experimental  values  from  Ref.  [54]. 


QCT  calculations 

Quasiclassical  trajectory  calculations  are  used  to  determine  thermal  rates  coefficients  from 

I - ^  ^ 

kiTt,  T„Tr,  Te)  =  J^  j  a{E,;  T*;  (14) 

’  0 

where  Tt,  T„,  and  Tg  are  the  translational,  vibrational,  rotational  and  electronic  tempera¬ 
ture,  jSi  =  l/{kBTi)^  i  denotes  the  respective  temperature,  Ec  is  the  collisional  energy,  and  ks 
is  the  Boltzmann  constant.  The  various  degrees  of  freedom  in  the  reactants  are  considered 
to  be  in  thermal  equilibrium,  i.e.  Tt  =  =  Tr  =  Te  and  it  will  be  referred  to  as  T  for  the 
remainder  of  the  work.  The  cross  section  a{Ec;  T)  (using  T  =  =  Tr  =  Te)  is  defined  by 


a{E,-,T) 


i^max  Jmax(j^) 

E  E 


u=0  j=0 


{2j  + 


Q 


rv 


(15) 


where  Qrv  is  the  rovibrational  partition  functions,  i'  and  j  are  the  vibrational  and  rotational 
quantum  numbers,  and  E^^j  is  the  respective  internal  energy  in  the  u  and  j  state.  The  state- 
to-state  cross  section  is 

CXD 

a{Ec;iy,j)  =  2tt  J  P^j{b;Et)bdb  (16) 

0 

where  Pyj{b;Et)  is  the  opacity  function,  and  b  the  impact  parameter. 


The  dynamics  of  the  system  is  followed  by  propagating  the  equations  of  motions  for  a  given  set 
of  initial  conditions.  These  equations  were  integrated  numerically  using  the  Velocity- Verlet 
algorithm  [55].  The  initial  conditions  for  NO  were  generated  from  a  WKB-quantized  periodic 
orbit  [56]  of  the  corresponding  rotating  Morse  oscillator  for  a  given  v  and  j  [57].  The  symmetry 
axis  of  NO,  the  axis  of  its  rotation,  and  the  angular  momentum  were  randomly  oriented  and 
the  integral  Eq.  14  was  evaluated  using  an  Importance  Sampling  Monte  Carlo  scheme.  [58] 
The  rotational  and  vibrational  quantum  numbers  and  the  collision  energy  Ec  were  sampled 
using  the  probability  distributions 


Puj  (T) 


(2j  -F  l)e(-^^-^') 


!^max  jmaxfj^O  _ 

E  E  (2/  +  l)e-^^-^V 


j'=0 


(17) 
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and 


piE,)dE,  =  (3^E,e-^^^dE,  (18) 

The  impact  parameter  was  sampled  from  0  to  6max  using  an  Importance  Sampling  Monte 
Carlo  schemefrom  which  the  rate  coefficient  is  determined  as 


k{T)  =  geiT) 


8  27r6ma; 
fi-rrfd  Ntot 


Nr, 


k=l 


(19) 


where  A^tot  and  A^reac  are  the  total  and  the  reactive  number  of  trajectories,  respectively,  and 
bk  is  the  impact  parameter  of  reactive  trajectory  k.  This  last  equation  includes  the  electronic 
degeneracy  factor  ge,  which  is  defined  as 

9e{T)  =  gN2o(^)%f4g)(r)g"Q(2n)(^)  (20) 

where  the  electronic  partition  functions  are  (?n20  =  3  (for  both  states  l^A'  and  l^A''), 
9n(4s)  =  4  and  gNO(2n)(7')  =  2  +  2 exp(-177.1/r). 


Overall,  10000  trajectories  at  nine  temperatures  between  100  K  and  20000  K  were  computed. 
For  each  temperature,  conservation  of  the  total  energy  was  verified.  A  time  step  of  0.01  fs  was 
used,  and  the  maximum  impact  parameter  was  set  large  enough  that  only  elastic  collisions 
take  place  (6m ax  =  10.05  A  and  9.58  A  for  the  ^A'  and  ^A"  respectively).  For  several  T,  a 
second  set  of  calculations  was  run  with  a  small  value  of  6niax  for  having  a  sufficient  number  of 
trajectories  to  analyse  the  exchange  channel.  The  final  rate  coefficients  were  calculated  from 
the  values  determined  in  the  calculation  on  each  PES, 

k{T)  =  kp{T)  +  kpp{T)  (21) 

where  kp  and  kpp  are  the  rate  coefficients  from  simulations  on  the  ^A'  and  ^A"  PES,  respec¬ 
tively. 


3.2  Results 

Quality  of  the  Potential  Energy  Surfaces 

The  quality  of  the  RKHS  PESs  is  reported  in  Figure  7.  The  cusp  regions  due  to  the  surface 
crossing  regions,  which  were  one  of  the  main  causes  of  errors  in  earlier  PESs  [47],  are  treated 
satisfactorily.  The  pronounced  angular  dependence  of  the  energy  as  a  function  of  R,  also  is 
shown  in  this  figure  for  a  typical  cut  at  tno  =  1-3  A. 

The  experimental  relative  energies  of  the  different  channels  can  be  compared  with  the 
analytical  PESs.  These  data  are  summarized  in  Table  2.  The  surfaces  developed  in  this  work 
underestimate  the  relative  experimentally  reaction  energy  (NO-I-N— 7-N2-I-0)  by  less  than  3.6% 
while  for  the  full  atomization  this  value  reduces  up  to  1.5%.  Li  et  al  [26]  used  an  extrapo¬ 
lation  procedure  for  matching  the  energies  of  the  reactant  and  products  to  the  experimental 
ones  which  is  empirical  and  probably  not  valid  globally.  Therefore,  we  decided  not  to  include 
any  semi-empirical  parameter.  The  approach  used  in  this  study  is  geared  towards  probing 
the  performance  of  a  given  level  of  theory  which  can  eventually  be  improved  by  morphing 
procedures  [59]  together  with  experimental  data  (e.g.  from  spectroscopy;  reactive  scattering. 
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Figure  7:  One-dimensional  potential  energy  curves  for  the  NO-I-N  grid  at  tno  =  1-3  A  for 
several  values  of  9  and  the  (left)  and  (right)  states.  The  lines  are  the  RKHS-energies 
and  the  points  correspond  to  the  ab-initio  energies. 


Table  2:  Reaction  energies  from  the  PESs  developed  in  this  work  and  the  corresponding 
experimental  data  (in  kcal/mol). 


reaction 

^A' 

^A" 

exp“ 

NO-tN^N-tN-tO 

154.89 

154.92 

152.60 

N2-tO^N-tN-tO 

228.02 

228.01 

228.40 

N2-tO^NO-tN 

73.12 

73.09 

75.80 

“  Experimental  values  from  [47]  taken  from  [60] 


1  2  3  1  2  3 

r^n  (A)  f fjo  (A) 


Figure  8;  Contour  plots  at  the  NNO  angle  of  110°.  Contour  are  incremented  in  step  of 
10  kcal/mol  (blue  lines  are  used  for  E  <  100  kcal/mol  and  red  line  for  E  >  100  kcal/mol). 
The  zero  of  energy  correspond  to  the  valley  of  the  N2-I-O  asymptotic  PES. 


etc.). 

The  NO-I-N  — )•  N2-I-O  is  a  direct  reaction  only  on  the  the  ^A"  surface  as  can  be  seen  in 
Figure  8.  Other  contour-plot  at  r^o  =  1-15  A  is  represented  in  Figure  9  for  the  ^A"  state. 
These  plots  also  show  the  smoothness  typical  of  the  kernel  interpolations. 

The  ^A'  surface  is  more  complex  than  the  ^A"  one.  In  this  case,  almost  all  the  minima 
and  transition  states  have  an  N-0  distance  close  to  1.3  A.  Therefore,  the  contour  plot  at 
tno  =  1-3  A  is  reported  in  Figure  10.  Similar  topography  was  shown  in  Figure  9  of  Ref.  [47]. 
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Figure  9:  Potential  energy  surface  for  NO+N  at  a  fix  diatomic  distance  tno  =  1-15  A  for 
the  ^A"  state.  Isocontours  are  shown  in  increments  of  4  kcal/mol  (blue  lines  are  used  for 
E  <  100  kcal/mol  and  red  line  for  E  >  100  kcal/mol).  The  zero  of  energy  correspond  to 
the  valley  of  the  N2+O  asymptotic  PES.  In  this  graphics,  the  linear  NON  corresponds  to  a 
Jacobi  angle  of  zero  degree 

This  agreement  with  Ref.  [47]  was  expected  as  the  ab-initio  calculations  in  both  studies  were 
done  at  the  MRCI  level,  even  though  in  the  present  study  the  energies  was  computed  using  a 
cc-pvtz  basis  set  and  represented  with  a  RKHS  method  while  Lin  et  al  employed  the  maug-cc- 
pVTZ  and  a  fitting  function  formed  by  a  sum  of  pairwise  term  and  permutationally  invariant 
polynomials  in  bond  orders.  The  main  difference  is  that  the  PESs  of  the  present  work  can  be 
employed  to  study  both,  NO+N  collisions  (including  atom  exchange)  and  N2  formation  for 
the  N2+O  channel.  The  long  range  part  of  the  potential  is  described  by  the  q^’^{Rk,  R)  kernel 
which  gives  the  correct  behavior  of  the  dispersion  and  induction  interactions  {R~^)  for  large  R. 


R  cos  0  (A) 


Figure  10:  Potential  energy  surface  for  NO+N  at  the  state  with  fixes  at  1.3  A.  At  this 
diatomic  distance  almost  all  stationary  points  can  be  seen.  Contour  are  incremented  in  step 
of  4  kcal/mol  (blue  lines  are  used  for  E  <  100  kcal/mol  and  red  line  for  E  >  100  kcal/mol). 
The  zero  of  energy  correspond  to  the  valley  of  the  N2+O  asymptotic  PES.  In  this  graphics, 
the  linear  NON  corresponds  to  a  Jacobi  angle  of  zero  degree. 


Dynamics 

The  PESs  were  used  to  follow  the  different  reaction  channels  over  temperatures  ranging  from 
100  K  to  20000  K.  Eigure  11  shows  the  distribution  of  probabilities  of  the  rovibrational  en¬ 
ergy,  vibrational  and  rotational  quantum  number  at  5000  K  for  both  products,  NO  and  N2, 
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separately.  It  can  be  seen  that  formation  of  N2  with  low  rovibrational  energies  is  allowed  on 
the  state  but  not  for  at  5000  K.  Consequently,  small  vibrational  quantum  numbers 
are  less  probable  in  the  ^A'  PES.  The  ^A"  surface  allows  direct  reactions  and  N2  can  finish  in 
the  ground  vibrational  state.  N2  formation  on  the  ^A'  PES  involves  several  stationary  points 
of  the  surface,  allowing  more  efficient  redistribution  of  energy  in  the  very  small  number  of 
trajectories  (32)  that  lead  to  N2.  Alternatively,  the  atom  exchange  reaction  (N^O+N^  — )■ 
N^O+N^)  on  the  ^A'  and  ^A"  surface  is  improbable  and  occurs  only  for  0.3  %  and  0.7  %  of 
the  trajectories  at  5000  K.  Therefore,  their  contribution  to  the  final  distributions  of  NO  is 
very  small  and  the  N^O+N^  — ?-N^O+N^  is  the  predominant  process  in  both  surface.  In  this 
case,  the  Ai^  =  0  are  the  most  favoured  transitions  with  a  probability  of  89%  in  the  ^A'  state 
and  82%  on  the  ^A"  surface.  Also  the  averaged  hnal  vibrational  energy  of  the  products  (NO 

and  N2),  defined  as  (Ei,)  =  are  represented. 
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Figure  11:  Rovibrational  distribution  at  T  =  5000  K. 

Figure  12  reports  the  probability  distributions  of  the  rovibrational  energies  together  with  the 
vibrational  and  rotational  quantum  number  for  simulations  at  T  =  20000  K.  Rovibrational 
states  with  energies  larger  than  the  dissociation  energy  of  the  diatomic  systems  were  found. 
At  this  temperature,  high  rotational  states  are  populated  and  the  effective  potential,  which 
includes  the  centrifugal  barrier  =  R(r)  +  ’  Prevents  dissociation.  Also,  full 

atomization  of  the  system  is  possible  but  only  a  few  trajectories  finished  in  this  channel.  As 
this  channel  is  not  of  particular  interest  to  this  work  it  was  not  further  considered.  The 
probability  of  the  N^O+N^  — ^  N^O+N^  channel  at  20000  K  is  3.1%  in  the  ^A'  state  and 
4.7%  for  the  ^A'b  These  values  show  a  pronounced  dependence  on  T  when  compared  with  the 
respective  probabilities  (<  0.7%)  at  5000  K.  The  barrier  in  the  surfaces  for  the  N^O+N^  — ^ 
N^O+N^  reaction  make  this  channel  significant  only  at  high  temperatures.  Figure  12  also 
shows  the  distribution  of  NO  due  to  the  N^O+N^  — ^  N^O+N^  and  N^O+N^  — ^  N^O+N^ 
channels.  If  no  rearrangement  occurs  during  the  collision,  the  rotational  distributions  have 
a  Boltzmann  behavior.  However,  the  formation  of  N2  and  the  exchange  of  nitrogen  do  not 
show  clearly  a  Boltzmann  distribution.  Therefore,  assuming  one  temperature  for  the  products 
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could  not  be  a  good  approximation. 
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Figure  12;  Rovibrational  distribution  at  T  =  20000  K.  In  parenthesis  the  number  of  trajec¬ 
tories  that  finish  in  each  channel. 

An  exploratory  set  of  time-independent  (TI)  quantum  reactive  scattering  calculations  were 
performed  using  both  PESs.  The  abc  code  [61]  was  used  for  studying  the  NO-I-N  collision 
at  several  total  energies  Et.  This  code  solves  the  time- independent  Schrodinger  equation  of 
a  triatomic  system  using  the  coupled-channel  hyperspherical  coordinate  method.  The  calcu¬ 
lations  were  performed  for  total  angular  momentum  J  =  0.  After  some  convergence  tests, 
the  maximum  value  of  j  included  in  any  channel  was  50,  and  the  maximum  internal  energy 
was  set  to  4.0  eV.  The  value  of  maximum  hyperradius  pmax  was  fixed  to  20  ao  and  400  step 
in  the  log-derivative  method  were  employed.  A  suitable  set  of  trajectories  including  all  pos¬ 
sible  states  at  each  of  the  four  analyzed  Et  was  computed  for  comparing  with  the  quantum 
calculations.  Figure  13  shows  the  vibrational  distribution  of  the  final  state  of  N2  using  TI 
and  QCT  methods.  The  large  discrepancies  were  found  at  Et  =  0.804  eV  using  the  ^A'  PES. 
However,  in  the  QCT  calculations  at  this  Et,  only  a  few  trajectories  were  reactive  and  the 
statistics  might  not  be  enough.  In  general,  the  agreement  is  quite  good. 

Table  3  presents  the  rate  coefficient  of  NO-I-N  computed  using  the  ^A''  and  ^A'  surfaces.  The 
rate  coefficients  for  the  formation  of  N2  increase  with  the  temperature  in  the  ^A'  state  faster 
than  in  the  ^A"  case.  That  is  due  to  the  presence  of  an  early  barrier  in  the  ^A'.  At  low 
temperatures  this  reaction  occurs  predominantly  on  the  ^A",  however,  at  higher  tempera¬ 
tures  (e.g.  10000  K),  the  contribution  of  both  PESs  is  comparable.  The  exchange  channel  is 
relevant  only  at  high  temperatures,  due  to  a  barrier  at  the  entry  of  this  channel  in  both  PESs. 


In  Figure  14,  the  total  rate  coefficients  for  the  NO-I-N — )-N2-|-0  reaction  with  others  values  in 
the  [1000;  5000]  K  interval  reported  in  the  literature  are  compared.  Furthermore,  the  experi¬ 
mental  values  recently  recommended  by  the  Jet  Propulsion  Laboratory  of  Pasadena  [62]  are 
included.  The  Gamallo  calculations  [45]  overestimate  the  experimental  values  from  1500  K. 
The  better  agreement  at  lower  temperatures  could  be  due  to  their  use  of  the  variational  tran- 
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Figure  13:  Vibrational  distribution  of  the  NO+N— 7-N2+0  reaction  using  time  indepen¬ 
dent  (TI)  and  QCT  methods  in  both  PESs,  ^A'  and  ^A". 


sition  state  theory.  At  low  temperatures,  the  inclusion  of  the  tunnelling  and  zero  points  of 
energies  effects  are  expected  to  be  relevant.  The  rates  coefficients  at  T  =  100  K,  200  K, 
and  300  K  were  computed,  but  the  discussion  in  this  work  will  turn  on  the  calculations  for 
T  >  1000  K.  Time  independent  studies  are  taking  place  at  the  moment  for  analysing  the 
influence  of  the  quantum  effects  at  low  temperatures  for  this  reaction.  However,  at  large  T 
the  quantum  effects  should  play  a  smaller  role,  and  the  rate  coefficients  can  be  determined 
from  QCT  calculations. 


Figure  14:  Rate  coefficient  of  the  N2  formation  at  low  temperatures.  Experimental  data  from 
Ref.  [62],  the  “Gamallo-2003”  values  were  taken  from  Ref.  [45]  and  “previous  PES”  was  also 
taken  from  Ref.  [45]  but  using  other  surfaces. 

The  behaviour  of  the  rate  coefficients  at  the  temperatures  relevant  for  the  hypersonic  flight 
regime  is  shown  in  Eigure  15.  The  expected  rate  increase  with  increasing  T  can  be  seen  in 
Figure  15.  Also,  the  exchange  channel  become  important  at  high  temperatures. 
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Table  3;  Rate  coefficients  in  cm^ molecules  ^s  ^  of  the  N(^S)+NO(^n)  reaction  at  several 
temperatures. 


r(K) 

fc(3A") 

fc(3A"+3A') 

N(4S)+N0(2n)^N2(iS)+0(3p) 
1000  9.16x10-14  1.08x10-11 

(1.09  ±  0.05)  xlO-ii 

1000“ 

2.32x10-13 

3.02x10-11 

3.05x10-11 

1000'’ 

8.10x10-13 

8.31x10-12 

8.32x10-12 

1000“ 

2000 

1.24x10-12 

1.64x10-11 

(2.32  ±  0.85)  X 10-11 
(1.76  ±  0.07)xl0-ii 

2000“ 

2.26x10-12 

3.07x10-11 

3.29x10-11 

2000'’ 

4.17x10-11' 

1.56x10-11 

1.60x10-11 

2000“ 

5000 

9.72x10-12 

2.96x10-11 

(2.24  ±  0.69)xl0-ii 
(3.93  ±  0.13)xl0-ii 

5000“ 

7.36x10-12 

4.06x10-11 

5.46x10-11 

5000'’ 

1.41x10-11 

3.72x10-11 

4.46x10-11 

10000 

3.13x10-11 

4.88x10-11 

(8.01  ±  0.23)  X 10-11 

15000 

4.36x10-11 

5.75x10-11 

(1.01  ±  0.03)xl0-i° 

20000 

5.46x10-11 

6.81x10-11 

(1.23  ±  0.03)xl0-i° 

N(4S)+N’0(2n)^N’(4S)+N0(2n) 
5000  7.22x10-13  1.41x10-12 

(2.13  ±  0.23)x10-i2 

5000“ 

1.19x10-12 

5.02x10-13 

1.69x10-12 

10000 

7.04x10-12 

1.14x10-11 

(1.85  ±  0.10)xl0-ii 

15000 

1.68x10-11 

2.19x10-11 

(3.87  ±  0.17)xl0-ii 

20000 

2.20x10-11 

2.75x10-11 

(4.94  ±  0.21)xl0-ii 

“  Value  computed  by  Gamallo  et  al  [45] 

^  Value  reported  in  Ref.  [45],  calculated  using  previous  PESs. 
°  Experimental  recommended  data  from  Ref.  [62] 
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T(K) 

Figure  15:  Rate  coefficient  for  the  NO+N  reaction  from  5000  K  to  20000  K.  The  points  are 
the  computed  rates  and  the  lines  corresponded  to  the  expressions  22  and  23. 


For  modelling  re-entry  conditions  the  use  of  an  expression  that  describes  the  rate  coefficients 
could  be  very  useful.  If  the  values  are  fitted  to  the  Arrhenius  expression  for  the  exchange  rate 
coefficients,  from  5000  to  20000  K  it  can  be 

J.exchfj,^  =  1.475  X  10“^°exp(-20907.67/r)  (22) 

In  the  case  of  the  N2  formation,  this  expression  from  5000  to  20000  K  have  the  form 

k{T)  =  1.712  X  10“^°exp(-7423.430/r)  (23) 

Furthermore,  a  more  general  expression  valid  in  a  largest  range  of  temperatures  (1000  K  to 
20000  K)  using  the  modified  Arrhenius  expression  can  be  written  as 

kmod{T)  =  2.421  X  exp(-102.447/r)  (24) 


3.3  Summary 

The  two  lowest  triplet  states  of  N2O  from  ab-initio  calculations  at  the  MRCI-I-Q  level  were 
presented.  The  energies  in  each  asymptotic  channel  were  fitted  using  a  kernel  based  method. 
Furthermore,  quasiclassical  trajectories  calculation  using  these  new  PESs  were  performed  up 
T  =  20000  K.  The  rate  coefficients  computed  in  this  work  were  compared  with  previous  studies 
up  to  5000  K.  At  low  temperatures,  the  quantum  effects  should  play  a  major  role  and  the  rates 
are  not  well  described  from  QCT  calculations.  However,  from  1000  K  the  results  of  this  study 
are  in  good  agreement  with  the  experimental  recommended  values.  At  temperatures  lowers 
than  5000  K  the  reaction  take  place  mainly  in  on  the  ^A"  surface.  For  higher  temperatures, 
the  contribution  of  both  states  should  be  considered.  Finally,  the  Arrhenius  expression  for 
the  rate  coefficient  at  the  extreme  temperatures  relevant  to  the  hypersonic  flight  regime  for 
the  N2  formation  and  the  exchange  channel  are  reported. 
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4  Collision-Induced  Rotational  Excitation  in  N^(^E+,v=0)— 
Ar:  Comparison  of  Computations  and  Experiment 

4.1  Introduction 

The  collisional  dynamics  of  cations  with  Ar  atoms  is  studied  using  quasi-classical 

simulations.  N^-Ar  is  a  proxy  to  study  cooling  of  molecular  ions  and  interesting  in  its  own 
right  for  molecule-to-atom  charge  transfer  reactions.  An  accurate  potential  energy  surface 
(PES)  is  constructed  from  a  reproducing  kernel  Hilbert  space  (RKHS)  interpolation  based  on 
high-level  ab  initio  data.  The  global  PES  including  the  asymptotics  is  fully  treated  within 
the  realm  of  RKHS.  Erom  several  ten  thousand  trajectories,  the  final  state  distribution  of 
the  rotational  quantum  number  of  after  collision  with  Ar  is  determined.  Contrary  to  the 
interpretation  of  previous  experiments  which  indicate  that  up  to  98  %  of  collisions  are  elastic 
and  conserve  the  quantum  state,  the  present  simulations  find  a  considerably  larger  number 
of  inelastic  collisions  which  supports  more  recent  findings. 

Cooling  neutral  and  charged  atoms  and  molecules  is  an  essential  step  for  controlled  investi¬ 
gations  of  fundamental  processes  in  chemistry  and  physics.  These  methods  rely  in  one  way 
or  another  on  energy  transfer  between  the  species  to  be  cooled  and  their  environment.  Many 
techniques  to  achieve  cooling  to  the  millikelvin  regime  exist.  They  include,  among  others, 
sympathetic  cooling  of  neutral  atoms  through  collisions  with  laser-cooled  species  [63]  or  cool¬ 
ing  through  collisions  with  Ar  atoms  in  crossed  molecular  beams.  [64] 

An  ionic  system  which  is  of  considerable  interest  in  this  context  is  N^-Ar.  It  is  a  prototypi¬ 
cal  system  to  study  molecular  ions  [65,  66]  and  charge  transfer  reactions  for  different  internal 
states  of  the  diatomic.  [67]  Previous  experiments  showed  that  the  charge  transfer  from  the 
molecule  to  the  atom  occurs  only  if  is  vibrationally  excited  and  no  charge  transfer  is  pos¬ 
sible  for  =  0).[65]  Figure  16  shows  relevant  stages  for  the  N^(u',y)  -|-  Ar  — )■  N2{v"J") 

-1-  reaction. 

Another  process  which  takes  place  on  the  electronic  ground  state  potential  energy  surface 
(PES)  is  scattering  of  the  Ar  atom  from  the  cation  through  formation  of  the  collision  com¬ 
plex  [N2Ar]^  (right  hand  side  in  Figure  16).  Recent  experiments  of  sympathetically  cooled 
N2+  (u  =  0,  =  0,  J  =  0.5)  and  collisions  with  room  temperature  Ar  atoms  have  shown 

that  the  measured  rate  for  quadrupole  vibrational  excitation  of  in  Coulomb  crystals  is 
two  orders  of  magnitude  lower  than  expected.  [71]  It  was  speculated  that  this  is  due  to  de¬ 
population  of  the  initially  prepared  molecular  state  through  rotationally  inelastic  collisions 
between  the  ion  and  the  Ar  buffer  gas.  This  process  has  been  previously  investigated 
using  laser  excited  in  a  22-pole  ion  trap  in  the  presence  of  a  low  density  Ar  buffer  at  90 
K.  Specifically,  the  rotational  state  J  =  6.5  of  a  small  number  of  N2-I-  ions  was  depopulated 
via  laser  induced  charge  transfer  with  subsequent  observation  of  refilling  the  hole  via  inelastic 
collisions.  Analysis  of  the  experimental  data  using  an  advanced  kinetic  model  yields  a  very 
low  inelastic  collision  rate  (<  2%). [65]  In  order  to  address  these  conflicting  experimental  re¬ 
sults  and  to  provide  a  more  atomistically  resolved  picture,  it  is  desirable  to  investigate  this 
system  using  accurate  computations.  Such  an  approach  has  provided  further  insight  into  re¬ 
lated  elementary  processes  including,  e.g.,  the  rotational  excitation  of  through  collisions 
with  N2,[72]  and  the  reactive  collisions  of  NO-|-O[10,  73]  and  OH-|-H[74]. 
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Figure  16:  Schematic  representation  of  the  relevant  stages  in  the  N2(u',/)  +  Ar"*"  — )•  {v'\j") 
+  Ar  reaction.  Experimentally  determined  values  are  indicated  by  “exp”  where  the  super¬ 
script  (a)  refers  to  Ref. [68]  and  (b)  to  Ref(s).[69,  70]  For  the  vibrational  and  rotational  ground 
state,  the  N^(u",y')  -|-  Ar  channel  is  0.179  eV  below  the  N2(u',y)  -t-  Ar+  asymptote. [69,  70] 
The  reported  experimental  dissociation  energy  of  the  [N2Ar]"*“  complex  ranges  from  1.04  to 
1.23  eV.[68] 


The  present  work  focusses  on  the  final  state  distribution  of  the  rotational  quantum  number 
j”  of  the  ion  after  collision  with  Ar  starting  from  a  given  initial  rotational  state  j'  (right 
hand  side  in  Figure  16).  Such  simulations  provide  information  about  the  conservation  of  the 
internal  state  of  during  and  after  collisions.  Because  of  the  appreciable  binding  energy 
of  1.14  eV  stabilizing  the  [N2Ar]"*'  complex  (see  Figure  16),  it  is  expected  that  the  dynamics 
in  the  bound  state  [N2Ar]'’'  leads  to  inelastic  collisions  and  gives  rise  to  rotational  excitation. 
Earlier  experimental  results  suggest  the  contrary,  namely  that  rotational  transitions  in 
-|-  Ar  collisions  occur  rarely  and  up  to  98  %  are  elastic  (“..implies  that  only  one  out  of  50 
collisions  results  in  a  change  of  the  rotational  state” ) .  [65]  A  possible  explanation  for  the  sur¬ 
prising  conservation  of  the  initial  rotational  state  could  be  the  consequence  of  some  hidden 
constants  of  motion  leading  to  approximate  selection  rules. [65] 

The  collision  of  molecules  with  Ar  atoms  is  studied  through  quasi-classical  molecular 
dynamics  (MD)  simulations.  Hence,  quantum  effects  due  to  zero-point  vibrations  and  tun¬ 
neling,  which  may  influence  the  collisional  dynamics  and  the  resulting  rotational  excitation 
of  N+,  are  not  included  in  the  simulations.  Moreover,  rotational  fine  structure  originating 
from  coupling  of  the  spin  of  and  the  angular  momentum  of  the  diatomic  are  not  included 
in  the  present  approach.  Technically,  accurate  full  close-coupling  quantum  scattering  calcu¬ 
lations  even  for  triatomic  systems  with  a  deep  well  (here  in  excess  of  1  eV),  large  anisotropy 
of  the  interaction  potential  and  small  rotational  constant  of  the  diatomic  are  currently  not 
feasible.  On  the  other  hand,  quasi  classical  trajectory  simulations  have  provided  valuable 
insight  into  rotational  excitation  in  the  related  charge  exchange  reaction  -|-  N2  — )•  N2  -|- 
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Figure  17:  +  Ar  system  described  in  Jacobi  coordinates.  The  dark  blue  spheres  A  and  B 

correspond  to  the  nitrogen  atoms  and  the  light  blue  sphere  C  to  the  argon  atom. 

N^.[72]  Similar  classical  studies  have  been  performed  for  the  O+CHD3  system[75]  and  the 
hydrogen  exchange  reaction  OH“+HBr  — )•  Br“+H20  without  accounting  for  spin-orbit  cou¬ 
pling  (as  in  the  present  case)  and  still,  good  agreement  with  experiment  was  obtainted.[76] 
In  the  light  of  the  computational  demands  of  a  fully  quantum  mechanical  treatment  and  the 
aforementioned  insights  obtained  for  similar  systems,  a  classical  treatment  appears  to  be  a 
meaningful  approach.  To  further  asses  the  role  of  quantum  effects,  the  classical  results  are 
complemented  with  calculations  accounting  for  quantum  zero  point  energy  (ZPE)  using  two 
different  methods. 

The  ZPE  motion  of  may  have  non-negligible  effects  on  the  results.  In  particular,  since  all 
simulations  are  carried  out  within  the  framework  of  classical  mechanics,  ZPE  might  be  con¬ 
sumed  and  lead  to  excitation  of  internal  degrees  of  freedom,  which  is  unphysical.  Since  there 
is  no  unique  way  to  correct  for  ZPE, [77]  two  different  strategies  to  probe  this  are  explored: 
a)  trajectories  are  run  on  the  bare,  ab  initio  PES  and  a  binning  criterion  is  applied  at  the 
analysis  stage,  and  b)  the  trajectories  are  run  on  a  ZPE-corrected  PES  in  order  to  prevent  the 
unphysical  consumption  of  ZPE  during  the  simulations.  Both,  a  binning  criterion[78]  and  a 
way  to  constrain  vibrational  motion[79,  80]  have  been  previously  employed  in  the  literature. 
The  PES  is  constructed  from  ab  initio  points  using  interpolation  by  the  reproducing  kernel 
Hilbert  space  (RKHS)  formalism.  [81,  82]  Since  the  present  work  focuses  on  the  study  of  ro¬ 
tational  transitions  in  (u^=0,/)  -|-  Ar  collisions,  only  the  electronic  ground  state  of  is 
considered. 


4.2  Methods 

Ab  initio  calculations 

The  N^-Ar  system  is  described  in  Jacobi  coordinates,  where  r  is  the  N-N  distance,  R  is  the 
distance  between  Ar  and  the  center  of  mass  of  the  nitrogen  atoms,  and  a  is  the  angle  between 
the  direction  of  R  and  the  N-N  axis  (see  Eigure  17). 

A  total  of  3025  points  on  a  regular  grid  was  used  to  sample  the  PES  at  the  UCCSD(T)  level 
of  theory  with  an  aug-cc-pVTZ  basis  set  [83,  84]  using  the  Gaussian  09  suite  of  codes.  [85]  The 
angular  grid  corresponds  to  an  11-point  Gauss-Legendre  quadrature  (a  [°]  =  11.815,  27.452, 
43.089,  58.726,  74.363,  90.000,  105.637,  121.274,  136.911,  152,548,  168.185). [86]  Eor  r,  the 
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grid  consists  of  11  points  according  to  the  equilibrium  position  req  of  the  bond  and  the 
turning  points[87]  from  u  =  0tou  =  4(r  =  0.998,  1.010,  1.024,  1.043,  1.071,  1.122,  1.167, 
1.207,  1.237,  1.263,  1.287  A).  The  grid  of  van-der-Waals  distances  R  was  chosen  to  capture 
the  well-region  with  more  densely  spaced  points  and  the  asymptotic  regions  with  fewer  points, 
for  a  total  of  25  grid  points  {R  =  1.7,  1.8,  1.9,  2.0,  2.1,  2.2,  2.3,  2.4,  2.5,  2.6,  2.7,  2.8,  2.9,  3.0, 
3.1,  3.2,  3.3,  3.4,  3.5,  3.6,  3.8,  4.1,  4.5,  6.0,  9.0  A). 

The  geometry  of  the  [N2Ar]'’'  complex  was  optimized  and  the  most  stable  configuration  was 
found  to  be  linear  with  Jacobi  coordinates  of  the  optimized  structure  of  req  =  1.105  A, 
i7eq  =  2.749  A  and  cteq  =  0°.  The  dissociation  energy  D^.  of  the  complex  in  the  linear  config¬ 
uration  is  1.266  eV  with  a  zero  point  energy  of  0.129  eV  (calculated  at  the  UCCSD(T)/aug- 
cc-pVTZ  level  of  theory).  Thus  the  dissociation  energy  Dq  of  1.137  eV  is  in  good  agreement 
with  the  experimentally  measured  value  of  1.109  eV.[88] 


RKHS  interpolation 

In  the  following  only  a  description  of  the  RKHS  formalism  is  provided  for  completeness.  For 
more  details  the  reader  is  referred  to  the  literature.  [89] 

General  procedure 

For  constructing  a  multi-dimensional  representation  of  a  PES  K(x)  for  a  molecular  system 
based  on  N  ab  initio  data  points  at  an  arbitrary  configuration  x,  V (x)  is  considered  to 
be  a  bounded  function.  The  quality  of  the  representation  of  K(x)  in  terms  of  a  kernel  is 
determined  by  two  aspects;  one  is  the  number  of  grid  points  which  equals  the  number  of 
kernel  coefficients  that  are  used,  and  the  other  is  the  type  of  kernel  function  itself.  As  such, 
V  (x)  can  be  represented  in  the  following  fashion 

N 

l/(x)  =  ^  WiQ{:>Ci,  x)  (25) 

i=l 

where  x^  are  the  coordinates  of  the  system  at  the  i-th  grid  point,  Wi  are  the  kernel  coefficients, 
and  (5(xj,x)  is  a  reproducing  kernel. 

The  multidimensional  reproducing  kernel  Q{x,  x')  can  be  represented  as  a  product  over  1- 
dimensional  kernels. [89]  For  example,  if  x  corresponds  to  the  Jacobi  coordinates  (r,  R,  a), 
the  3-dimensional  kernel  can  be  expressed  as 

(5(x,x')  =  qi{r,r')  ■  q2{R,R')  ■  q3ia,a).  (26) 

where  the  kernels  qi  can  be  chosen  freely  depending  on  the  nature  of  the  respective  variable 
(angular  or  distance-like).  The  kernel  coefficients  Wi  in  Eq.  25  are  determined  by  solving  the 
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following  system  of  linear  equations, 


Q(xi,Xi) 

Q(xi,X2)  . 

.  .  Q(xi,XAr) 

Wl 
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(5(x2,X2)  . 
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(5(xAr,X2)  . 
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_WN_ 

"  initio 

(xi)‘ 

initio 

(X2) 

V 


ab  initio 


(XTV). 


(27) 


Since  the  reproducing  kernel  matrix  is  symmetric  and  positive  definite  by  definition,  the  com¬ 
putationally  efficient  Cholesky  decomposition  can  be  used  to  solve  Eq.  28.  [90]  Once  the  kernel 
coefficients  Wi  are  known,  the  energy  at  an  off-grid  point  x  can  be  evaluated  from  Eq.  25. 


1-Dimensional  Kernels 

Distance-like  coordinates:  In  this  work,  the  1-dimensional  kernel  for  distance-like  coordinates 
was  that  from  Ref.  [89], 

(28) 

xB(m-\-  1,  n)2Ei(— n  — ) 

x> 

where  B{m  -\-  1,  n)  is  the  beta  function,  2Fi{—n  -|-  1,  m  -|-  1;  n  -|-  m  -|-  1;  x</x>)  is  the  Gauss 
hypergeometric  function  and  x<  and  x>  denote  the  smaller  and  the  larger  of  x  and  x'  respec¬ 
tively.  The  chosen  value  for  n  controls  up  to  which  derivative  the  kernel  is  smooth,  while  the 
value  of  m  controls  its  long-range  behaviour.  [89]  In  particular,  if  no  points  are  available  in 
the  asymptotic  region,  m  can  be  chosen  to  mimic  the  physical  long-range  behaviour  of  the 
interpolated  variable.  Here  we  use  n  =  2  and  m  =  6  for  both  distance-like  kernels,  although 
different  values  are  possible.  [91]  In  the  present  work  the  choice  of  m  is  largely  inconsequential 
because  enough  data  is  available  in  asymptotic  regions  and  the  ab  initio  points  are  necessarily 
reproduced  exactly. 

Angle-like  coordinates:  Eor  angle-like  internal  coordinates  the  general  expression  for  the  kernel 
is  [89] 

n—  1 

q^{x,x')  =  ^x>x<  -I- nx<x>“^2Ei(l, -n -I-  l;n-|-  1;  — )  (29) 

i=o 

Eq.  29  is  valid  only  if  the  angle-like  coordinate  is  rescaled  so  that  both  x  and  x'  belong  to  the 
interval  [0, 1].  Eor  example,  to  rescale  the  Jacobi  coordinate  a  a  new  coordinate  y  is  defined 
as  y  =  (1  —  cosq;)/2.  Again,  choosing  n  =  2  is  sufficient  for  an  accurate  representation  of  the 
angular  dependence. 

Treatment  of  the  Asymptotics 

The  1-dimensional  distance-like  kernel  given  by  Eq.  29  has  the  advantage  that  it  correctly 
decays  to  zero  at  long  range.  While  this  is  still  true  in  the  present  case  for  one  particular  value 
of  the  N-N  separation  r,  this  is  not  the  case  for  the  total  energy  of  the  complex  for  arbitrary 
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Figure  18:  Contour  plot  of  the  RKHS  interpolated  N^-Ar  PES  (white  lines)  with  at  its 
equilibrium  position  (igg  =  1.122  A).  The  binding  energy  is  around  1.3  eV  and  isocontours 
are  labelled  with  energies  in  eV.  The  color  map  is  the  error  surface  with  relative  errors  in 
%.  Relative  errors  were  obtained  by  comparing  the  interpolated  energies  to  off-grid  ab  initio 
data.  Even  the  largest  relative  errors  are  below  0.1  %. 


values  of  r.  The  correction  of  the  asymptotics  was  considered  in  previous  work  on  HO2  by 
using  a  manybody  expansion  of  the  PES.  [92]  For  the  present  work,  a  different  method,  which 
is  fully  within  the  RKHS  formalism,  was  used  and  is  presented  in  Section  4.3. 

Correction  of  ZPE  in  the  AJ -Ar  PES  One  way  to  account  for  ZPE  in  the  simulations  is 
to  carry  out  a  point-wise  correction  of  the  PES  whereby  the  ZPE  contribution  is  taken  into 
account  at  each  configuration.  [93]  This  contribution  is  calculated  from  the  harmonic  oscillator 
energy  evaluated  for  n  =  0 

yZPE  (30) 

i=l  ^  ^ 

where  A  =  3  is  the  number  of  degrees  of  freedom  of  the  system  and  Wj  are  the  positive 
eigenvalues  obtained  from  ab  initio  calculations  at  the  MP2/aug-cc-pVTZ  level  of  theory. 

In  the  present  simulations  two  different  PESs  were  employed,  one  where  the  ZPE  correction 
given  by  Eq.  30  is  included  [ZPE- corrected  PES)  and  one  where  no  such  correction  was  made 
[bare  PES).  The  RKHS  PES  and  its  quality  compared  to  off-grid  points  is  shown  in  Figure 
18. 
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Molecular  Dynamics  Simulations 

MD  simulations  of  the  collision  between  and  Ar  have  been  carried  out  with  CHARMM.[94] 
The  equations  of  motion  were  integrated  nsing  a  time  step  of  At  =  0.1  fs  and  energies  were 
calculated  from  the  RKHS  representation  of  the  PES  described  in  Section  4.2  using  the  bare 
and  ZPE-corrected  PESs,  respectively.  Forces  were  obtained  by  numerical  differentiation  us¬ 
ing  a  five-point  stencil.  [95] 


time  [ps] 


Figure  19:  A  representative  molecnlar  dynamics  simulation  for  the  collision  between  and 
Ar.  The  upper  panel  reports  the  variation  of  the  total  energy  as  a  function  of  time  for  time 
steps  At  =  0.5  and  0.1  fs,  respectively.  In  the  bottom  panel  the  distance  between  the  ion  and 
the  incoming  Ar  atom  is  shown  for  At  =  0.1  fs.  After  ~  0.5  ps  the  complex  is  formed  and 
lives  for  about  12.5  ps  after  which  it  breaks  up  again. 

Initial  conditions  were  generated  from  a  WKB-qnantized  periodic  orbit  of  the  corresponding 
rotating  Morse  oscillator  for  given  vibrational  v  and  rotational  j  quantnm  numbers.  [16,  96] 
Parameters  for  the  Morse  potential  were  fitted  to  ab  initio  data  for  {De  =  8.5003  eV, 
Te  =  1.1307  A,  P  =  2.5928  A“^).  The  fitted  dissociation  energy  is  consistent  with  the  exper¬ 
imental  value  of  8.724  eV.[97]  The  initial  vibrational  quantum  number  was  v'  =  0  and  the 
rotational  quantum  numbers  were  either  f  =  0  or  j'  =  6.  Relative  collision  energies  between 
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and  Ar  were  sampled  from  a  Maxwell-Boltzmann  distribution  at  either  90  K  or  300  K. 
These  initial  conditions  are  representative  of  the  previous  experiments  which  investigate  [N^ 
{v'  =  0,/  =  6);  T  =  90  K][65]  and  the  collision  of  ultracold  molecular  ions  with  Ar  at  room 
temperature  [N^  {v'  =  0,j'  =  0);  T  =  300  K].[71]  The  impact  parameter  b  was  selected  from 
a  uniform  distribution  between  0  and  25  oq  (prior  test  runs  had  shown  that  larger  impact 
parameters  do  not  lead  to  collision  or  even  mutual  influence  of  the  impact  partners).  Such 
a  procedure  was  already  employed  in  previous  studies  of  the  collisions  between  0(^P)  and 

NO(2n).[io] 

Simulations  on  the  ZPE-corrected  PES  used  a  suitably  modified  procedure  for  generating 
initial  conditions  for  the  positions  and  velocities.  Eor  u'  =  0  and  any  j' ,  the  position  is  always 
set  to  the  minimum  of  the  effective  Morse  potential  14fr(^)  =  V{r)  —  ^  and  all  vibrational 
energy  is  removed  from  the  N-N  stretch  as  the  ZPE  is  already  accounted  for  in  the  corrected 
PES.  The  vibrational  energy  corresponding  to  the  ZPE  is  added  back  prior  to  the  analysis  in 
order  to  be  able  to  use  the  same  filtering  criteria  for  the  ZPE-corrected  and  bare  PES. 

All  trajectories  were  run  until  the  collision  partners  were  fully  separated,  for  which  a  value 
corresponding  to  more  than  1.3-times  their  initial  separation  was  assumed.  The  maximum 
simulation  time  considered  was  1  ns.  For  statistically  significant  results  a  total  of  25000 
trajectories  was  run  for  each  temperature  and  PES  and  approximately  4000  of  them  were 
excluded  from  the  analysis  due  to  lifetimes  longer  than  1  ns. 


4.3  Results 

Correction  of  the  asymptote 

Depending  on  the  bond  length  r  individual  manifolds  V{R,a;r)  dissociate  to  different 
asymptotic  values  for  i?  — )■  oo.  This  assumes  that  the  PES  becomes  isotropic  (a— independent) 
with  increasing  R  which  is  in  general  true.  In  order  to  employ  RKHS  interpolation  in  a  mean¬ 
ingful  fashion  for  the  R—  and  a—  degrees  of  freedom,  the  asymptotic  value  for  each  manifold 
characterized  by  a  specific  value  for  the  N-N  separation  must  be  shifted  to  zero  energy.  This 
is  necessary  because  kernels  decay  to  zero  asymptotically. 

In  order  to  set  the  energy  to  zero  for  every  cut  (r  =  constant),  the  energy  of  the  isolated 
diatomic  E(r,R  — )•  oo)  is  subtracted  from  the  energy  of  all  points  that  share  the  same  value 
of  r  and  yields  i?,  a)  =  E{r,R,a)  —  E{r,R  — >  oo)  and  V°°{r)  =  E{r,R  — >  oo).  The 

two  data  sets  i?,  a)  and  V°°{r)  are  then  interpolated  within  the  RKHS  framework. 

Because  the  PES  is  isotropic  for  sufficiently  large  R,  the  interpolation  of  P°°(r)  requires  only 
a  I-dimensional  kernel.  Hence,  all  quantities  required  are  represented  by  kernels  and  the 
global  PES  can  be  evaluated  at  an  off-grid  point  (r^,  R^,  a^)  according  to  Eq.  31. 

P(r°,R°,a°)  =  P“'"(r°,R°,a°)  +  H“(r°)  (31) 
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Figure  20:  Part  of  the  global  PES  for  a  triatomic  system  N^-Ar.  The  relevant  coordinates 
(r,  R,  a)  are  indicated  in  the  inset  with  the  dark  blue  spheres  corresponding  to  the  molecule 
and  the  light  blue  sphere  to  the  Ar  atom.  The  three  curves  refer  to  different  intermolecular 
distances  of  the  diatomic  (^nn))  ^"nn  =  1-00  A  (red),  tnn  =  1-11  A  (blue)  and  tnn  = 
1.21  A  (black).  In  the  quantum  chemistry  calculations  the  asymptotes  of  each  manifold 
characterized  by  the  particular  tnn  separation  dissociates  to  its  own  asymptote  as  indicated 
by  the  non-overlapping  curves  at  long  range.  If  the  asymptote  is  not  corrected  in  the  RKHS 
interpolation,  all  curves  converge  to  a  value  of  E  =  0  for  large  R. 
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Simulation  and  Analysis  of  Collisional  Excitation 

Simulations  on  the  bare  and  ZPE-corrected  PES  were  carried  out  with  a  range  of  initial  con¬ 
ditions.  The  physical  conditions  included  two  temperatures,  T  =  90  K  and  T  =  300  K.  The 
first  corresponds  to  previous  experiments  [65]  and  the  latter  was  chosen  to  examine  the  role  of 
temperature  on  the  results.  The  time  step  in  all  simulations  was  At  =  0.1  fs.  Overall,  25000 
trajectories  were  run  for  each  temperature.  However,  a  considerable  number  of  trajectories 
are  not  further  analyzed  because  they  violate  in  one  way  or  another  restrictions  imposed 
by  quantum  mechanics.  Such  filtering  at  the  post-processing  stage  is  a  usual  procedure  to 
limit  the  number  of  classical  trajectories  to  be  analyzed  to  those  which  correspond  to  valid 
simulations  within  a  semiclassical  framework. 

Filtering  at  the  post-processing  stage  was  carried  out  as  follows.  The  WKB  procedure  can  be 
used  to  determine  the  vibrational  quantum  number  v"  on  the  effective  potential.  In  general  a 
real,  non-integer  number  is  obtained  for  v".  By  applying  a  binning  criterion,  only  trajectories 
with  v”  within  certain  thresholds  are  retained  for  further  analysis.  This  ensures  that  only 
trajectories  are  analyzed  in  which  ZPE  was  not  transformed  to  other  forms  of  energy  during 
the  simulation.  Three  different  thresholds  were  used  here:  frac(u'')  =  ±0.1,  ±0.01  and  ±0.001, 
where  frac(x)  is  the  difference  of  x  to  the  closest  integer  value.  These  values  correspond  to  a 
conservation  of  the  ZPE  within  220.70  cm“^,  22.07  cm“^  and  2.21  cm“^,  respectively. 

Rotational  j”  states  are  rounded  to  the  next  closest  integer  value  divisible  by  2  in  order  to 
fulfill  restrictions  imposed  by  the  symmetry  of  N^,  which  can  exist  in  either  the  ortho  or  para 
nuclear  spin  state.  If  a  binning  criterion  similar  to  the  vibrational  quantum  number  is  ap¬ 
plied  to  j”  states,  the  number  of  trajectories  that  fulfill  both  criteria  decreases  dramatically. 
Hence,  no  such  binning  was  applied  for  j” .  For  the  final  state  analysis  of  trajectories  run  on 
the  corrected  PES,  the  ZPE  is  added  back  into  the  classical  rovibrational  energy,  since  it  was 
removed  before  generating  the  initial  conditions. 

The  final  distribution  of  j"  states,  is  determined  from  the  WKB  procedure  after  the 

complex  has  dissociated.  An  event  is  classified  as  “inelastic”  if  the  final  j”  differs  from  the 
initial  j' ,  i.e.  j”  ^  j' .  As  the  simulations  involve  numerical  imprecisions  of  about  3.5  cm“^ 
(see  Figure  19,  upper  panel),  corresponding  to  roughly  the  energy  difference  between  /  =  0  to 
j"  =  1  for  (based  on  a  rotational  constant  oi  B  =  1.932  cm“^)[98],  an  additional  criterion 
on  the  necessary  amount  of  change  in  j”  was  introduced  to  classify  inelastic  transitions.  For 
this,  the  quantity  j*  is  used.  Only  trajectories  for  which  the  final  j"  differs  by  more  than  j* 
from  the  initial  j'  state  are  considered  to  correspond  to  an  inelastic  collision. 

Rotational  excitation  Figure  21  reports  the  final  state  distribution  P{j”)  for  from  sim¬ 
ulations  on  the  bare  and  ZPE-corrected  PES  starting  from  j'  =  0.  Results  are  presented  for 
simulations  at  T  =  90  K  and  T  =  300  K.  The  distributions  differ  slightly  between  the  two 
PESs  while  temperature  has  only  a  minor  effect  on  P{j").  On  the  bare  PES  higher  j"— states 
are  populated  compared  to  simulations  on  the  ZPE-corrected  PES.  Possibly,  this  is  due  to 
the  more  isotropic  shape  of  the  ZPE-corrected  PES  due  to  angular  averaging.  Also,  as  the 
well  is  more  shallow  on  the  ZPE-corrected  PES  due  to  inclusion  of  ZPE,  the  anisotropy  near 
the  inner  wall  is  less  accessible.  The  maximum  of  the  distribution  is  at  j"  =  2,  whereas  it  is 
at  j"  =  0  for  the  bare  PES.  This  already  suggests  that  more  than  a  fraction  (i.e.  2%  in  the 
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experiments [65])  of  the  collisions  are  inelastic. 


The  distribution  changes  only  slightly  between  different  filtering  criteria,  as  is  shown  in  the 
inset  in  Figure  21.  The  percentage  of  inelastic  and  elastic  collisions  for  the  bare  PES  and 
a  filtering  criterion  of  frac(u'')  =  ±0.01  is  summarized  in  Table  4.  A  collision  is  inelastic  if 
the  rotational  quantum  number  of  the  diatomic  differs  from  the  initial  j'  by  more  than  j* , 
i.e.  j”  ^  [/  —  j* ,j'  ±  j*].  Even  with  an  unrealistically  large  margin  of  j*  =  4,  the  fraction  of 
inelastic  collisions  is  well  above  the  experimentally  reported  value  of  2  %.[65]  Depending  on 
the  figure  of  merit  j*  used,  20  %  to  40  %  of  collisions  are  inelastic,  see  Table  4. 

Figure  22  shows  the  correlation  between  impact  parameter  b  and  j”  states.  It  is  typically 
found  that  smaller  impact  parameters  lead  to  higher  rotational  excitation  whereas  collisions 
with  large  b  are  preferentially  elastic. 

A  remarkable  result  is  that  excitation  to  the  highest  rotational  states  is  only  found  for  trajecto¬ 
ries  with  shorter  lifetimes  (see  Figure  23).  Closer  examination  of  some  individual  trajectories 
with  short  lifetime  shows  that  no  tight  [N2Ar]'*“  complex  has  to  be  formed  for  rotational  exci¬ 
tation  of  to  occur.  A  “tight  complex”  refers  to  a  situation  in  which  the  collision  partners 
come  close  enough  to  at  least  once  enter  the  short-range  repulsive  region  of  the  interaction 
potential.  For  rotational  excitation  it  is  sufficient  for  the  two  collision  partners  to  fly  past 
each  other  such  as  to  influence  their  respective  flight  paths.  The  interaction  between 
and  Ar  leads  to  a  torque  which  results  in  pronounced  rotational  excitation  of  N^.  Such 
^^fly-by”  trajectories  (lifetime  of  the  complex  shorter  than  5  ps)  lead  to  excitation  to  higher 
rotational  states  than  trajectories  with  a  longer  lifetime  and  make  up  21.5  %  of  all  trajectories. 

Counter-intuitively,  the  impact  parameter  b  is  not  correlated  to  the  lifetime  of  the  complex 
(see  inset  of  Figure  22).  In  fact,  the  complex  can  be  formed  and  live  for  a  very  long  time  even 
for  large  impact  parameters,  provided  that  the  collision  energy  is  sufficiently  low.  Conversely, 
small  impact  parameters  still  sometimes  lead  to  fly-by  trajectories,  provided  that  the  collision 
energy  is  sufficiently  large.  It  is  the  combination  of  a  small  impact  parameter  and  a  high 
collision  energy  (leading  to  fly-by  trajectories)  that  leads  to  the  highest  rotational  excitations. 


elastic  (%) 

inelastic  (%) 

T  =  90  K  ij*=2) 

bare  PES 

63 

37 

T  =  300  K  (V=2) 

67 

33 

T  =  90  K  (i*=4) 

70 

30 

T  =  300  K  (V=4) 

79 

21 

T  =  90  K  {j*=2) 

ZPE-corrected  PES 

64 

36 

T  =  300  K  (V=2) 

60 

40 

T  =  90  K  (i*=4) 

86 

14 

T  =  300  K  (V=4) 

78 

22 

Table  4;  Percentage  of  elastic  and  inelastic  collisions  for  a  hltering  criterion  of  frac(u')  =  0.01 
with  initial  j  =  0.  A  collision  is  considered  as  inelastic  when  j'  changes  from  0  to  j"  >  j* . 
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Figure  21:  Green  diamonds  (ZPE-corrected,  300  K),  blue  triangles  (ZPE- corrected,  90  K), 
black  circles  (bare,  300  K),  red  squares  (bare,  90  K).  Distribution  of  j"  states  after  dissociation 
for  ZPE  conservation  criterion  frac(u^')  =  ±0.01  and  different  temperatures.  Trajectories  were 
started  with  an  initial  j'  =  0  and  evolved  on  the  bare  and  ZPE-corrected  PES.  The  total  count 
of  trajectories  that  meet  the  conservation  criterion  are  4207  (T  =  300  K)  and  3309  (T  =  90 
K)  on  the  corrected  and  599  (T  =  300  K)  and  422  (T  =  90  K)  on  the  bare  PES,  respectively. 
Note  that  trajectories  on  the  corrected  PES  are  much  more  likely  to  meet  the  criterion, 
because  ZPE  is  removed  during  the  dynamics  and  added  back  in  the  analysis.  As  such, 
the  conservation  criterion  on  the  corrected  PES  can  only  be  violated  by  an  uptake  of  energy, 
whereas  it  can  be  violated  by  an  uptake  or  loss  of  energy  on  the  bare  PES.  Inset:  Distribution 
of  j"  states  after  dissociation  for  different  ZPE  conservation  criteria  (frac(u'')  =  ±0.1,  black 
circles,  frac(u")  =  ±0.01,  red  diamonds,  frac(u'')  =  ±0.001,  green  squares)  at  T  =  300  K  on 
the  ZPE-corrected  PES.  The  number  of  trajectories  that  meet  the  criterion  are  19821,  4207 
and  745,  respectively.  The  distribution  changes  slightly,  but  shows  the  same  overall  behaviour 
for  the  different  filtering  criteria. 
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Figure  22:  Impact  parameter  b  vs.  j"  {j'  =  0,  T  =  300  K)  on  the  corrected  PES.  Large  rota¬ 
tional  excitation  is  observed  only  for  small  impact  parameters.  Plots  for  other  temperatures 
on  the  corrected  and  bare  PES  show  similar  behaviour.  Impact  parameter  b  vs.  lifetime.  The 
lifetime  is  not  correlated  to  the  impact  parameter  (inset). 


lifetime  vs.  7 

3(X)  K.  ZPE-corrected) 


Figure  23:  Lifetime  vs.  j”  {j'  =  0,  T  =  300  K)  on  the  corrected  PES  for  different  ZPE 
conservation  criteria,  numbers  in  parentheses  indicate  how  many  trajectories  meet  the  crite¬ 
rion.  The  largest  rotational  excitation  is  observed  only  for  short  lifetimes.  Plots  for  other 
temperatures  on  the  corrected  and  bare  PES  show  similar  behaviour. 
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The  earlier  experiments [65]  were  thermal  and  started  from  a  distribution  of  rotational  states, 
but  the  initial  rotational  state  studied  was  j'  =  6.  Hence,  additional  simulations  were  carried 
out  for  this  case.  25000  independent  trajectories  were  started  on  the  ZPE-corrected  PES. 
Figure  24  shows  the  distribution  of  j"  states  after  the  complex  dissociates.  The  percent¬ 
age  of  inelastic  collisions  for  the  strictest  filtering  criterion  frac(u'')  =  ±0.001  and  j*  =  0 
is  41  %  and  decreases  to  17  %  for  j*  =  2,  respectively.  The  rate  constant  for  rotational 
excitation/relaxation  was  calculated  according  to  [10] 


k{T) 


8  27rb 


-^reac 

'max  7 

>  bi 


TTfiP  g{T)Ntot 


(32) 


where  /U  is  the  reduced  mass  of  the  N^— Ar  complex,  g{T)  is  the  electronic  degeneracy  factor, 
femax  is  the  maximum  impact  parameter,  bi  the  impact  parameter  of  trajectory  i,  iVtot  the 
total  number  of  trajectories  that  meet  the  filtering  criterion,  Nj-eac  the  number  of  inelastic 
trajectories  {\j''  —j'\  >  j*)  and  j3  =  l/{kBT).  The  calculated  value  of  =  1.17-10“® 

cm^  s“^  is  two  orders  of  magnitude  larger  than  the  reported  value  of  k  =  (1.4±0.4)  •  10“^^cm^ 
s“^.[65]  However,  it  should  be  noted  that  the  experimental  value  was  not  directly  measured, 
but  inferred  from  a  kinetic  model  based  on  the  rate  coefficient  for  charge  transfer  and  the 
rate  of  laser  excitation. 


Since  the  Langevin  rate  is  often  considered  to  be  an  upper  bound  for  reaction  rates  it  is 
surprising  that  the  computationally  determined  rate  kji=Q^jiijLQ  is  about  60  %  larger  than  the 
Langevin  rate  ki  =  7.4  •  10“^®  cm^  s“^.[65]  However,  rates  larger  than  kL  have  been  reported 
previously  in  the  literature  [99].  It  should  be  noted  that  Langevin  theory  assumes  an  idealized 
form  for  the  centrifugally  corrected  interaction  potential  between  the  ion,  which  is  modelled 
as  a  point  charge,  and  the  neutral  atom  given  by[100] 

1  a'e^ 

2  SttCoR'^ 

where  L  =  fivb,  v  is  the  relative  collision  velocity,  a'  the  polarizibility  volume  of  the  neu¬ 
tral  species,  /x  the  reduced  mass,  b  the  impact  parameter  and  R  the  distance  between  centre 
of  mass  of  the  molecular  ion  and  neutral  species.  A  comparison  between  the  1/R^  term  in 
Eq.  33  and  the  actual  interaction  potential  shows  that  Langevin  theory  is  insufficient  to 
describe  the  collisional  rate.  In  particular,  the  true  PES  decays  more  slowly  to  zero  and 
the  anisotropy  of  the  PES  (see  Figure  18),  which  is  crucial  for  the  dynamics,  is  completely 
neglected  in  the  Langevin  model.  Nonetheless,  if  Langevin  theory  is  applied  naively  to  the 
same  set  of  trajectories  that  was  used  to  calculate  the  /cj/=6-5-/V6  Eq.  32  is  used 

to  calculate  a  rate  (counting  those  trajectories  as  ’’reactive”  that  satisfy  Ec  >  VL(i?max), 
where  Ec  is  the  collision  energy  and  Hr(7?max)  is  the  Langevin  interaction  potential  at  the 
position  of  the  maximum  of  the  centrifugal  barrier  according  to  Langevin  theory),  a  rate 
constant  of  ki  =  6.2  •  10“^®  cm^  s“^  is  obtained.  This  is  in  good  agreement  with  the  value 
of  kc  =  7.4  •  10“^®  cm^  s“^  from  Langevin  theory.  Note  however  that  this  analysis  can  be 
performed  from  the  initial  conditions  without  running  actual  dynamics  and  merely  shows  that 
the  initial  conditions  for  our  simulations  are  consistent  with  the  Langevin  model.  However, 
the  model  itself  is  insufficient  to  describe  the  actual  dynamics.  It  should  also  be  considered 
that  ki  measures  merely  a  sort  of  “collision  rate” ,  but  as  was  pointed  out  earlier,  the  complex 
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does  not  need  to  be  formed  in  order  for  rotational  relaxation  or  excitation  to  occur.  In  fact, 
for  92.4  %  of  the  trajectories  the  Langevin  model  correctly  predicts  whether  the  complex  is 
formed  or  not  from  just  the  initial  conditions.  Since  the  observed  rate  is  60  %  larger  than 
the  Langevin  rate,  this  further  indicates  that  rotational  excitation  can  occur  without  complex 
formation.  A  final  test  was  to  run  simulations  with  an  explicitly  isotropic  PES  beyond  8  A.  At 
this  distance,  the  PES  is  still  appreciably  anisotropic  (~  50  cm“^  between  linear  and  T-shape 
geometry).  Hence,  the  PES  was  multiplied  by  an  empirical  factor  of  exp(— (i2/7.3A)^^),  which 
ensures  a  smooth  cutoff  at  long-range,  yet  leaving  the  short-range  part  of  the  PES  largely 
unaffected.  Although  this  PES  fulfills  the  requirement  of  Langevin  theory,  namely  that  the 
long  range  part  of  the  PES  should  be  isotropic,  the  computed  rate  is  unaffected  and  still 
exceeds  the  Langevin  rate.  It  is  suspected  that  the  anisotropy  in  the  range  below  8  A  is  the 
main  cause  for  rotational  excitation. 


corresponding  rotational  eneigy  [cm  '[ 


0  11.59  38.M  81.13  139.09  212.50  301.35  405.67  525.44  660.66 


Figure  24:  Distribution  of  j"  states  after  dissociation  for  different  ZPE  conservation  criteria  at 
T  =  90  K.  The  numbers  in  brackets  show  how  many  trajectories  meet  the  criterion.  For  every 
value  of  j,  the  corresponding  rotational  energy  is  given  (the  rotational  constant  is  B  =  1.932 
cm“^). [98]  Trajectories  were  started  with  an  initial  j'  =  6  and  evolved  on  the  ZPE-corrected 
PES. 

4.4  Conclusion 

Classical  molecular  dynamics  simulations  of  the  nonreactive  collision  between  the  cation 
and  Ar  atoms  at  two  different  temperatures  show  that  inelastic  rotational  excitation  of  the 
ion  in  the  product  channel  is  important  and  occurs  more  frequently  than  previously  assumed. 
The  simulations  use  an  RKHS  PES  based  on  UCCSD(T)/aug-cc-pVTZ  electronic  structure 
calculations  and  correct  handling  of  the  asymptotics  within  the  RKHS  framework.  Analysis  of 
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the  results  for  f  =  6  using  a  strict  filtering  criterion  of  frac(u'')  =  ±0.001  and  a  figure-of-merit 
j*  =  2  suggests  that  inelastic  collisions  occur  in  at  least  17  %  of  the  cases  which  is  one  order 
of  magnitude  larger  than  reported  in  earlier  experiments  (2  %).[65]  Interestingly,  the  [N2Ar]'’' 
complex  does  not  need  to  be  formed  (and  stabilized)  for  rotational  excitation  to  occur.  A 
sufficiently  close  encounter  of  the  two  collision  partners  is  sufficient  to  mutually  influence 
their  flight  paths  and  lead  to  rotational  excitation.  It  should  be  pointed  out  that  the  PES 
used  in  this  work  was  calculated  using  a  single-reference  method.  Electronic  effects,  which  are 
not  adequately  captured  using  single-reference  methods,  might  play  a  non- negligible  role  in 
the  dissociative  region  of  the  PES.  Consequently,  further  investigations  should  employ  multi¬ 
reference  methods  such  as  MRCI  to  capture  electronic  effects  which  are,  however,  outside  the 
scope  of  the  present  work.  For  a  complete  understanding  of  the  rate  of  rotational  excitation 
in  the  N^-Ar  system,  new  experiments,  which  allow  precise  control  of  the  exact  quantum 
state  of  the  collision  partners  and  additional  computational  investigations  at  the  quantum 
level  are  necessary.  [71] 
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